## bicameral conflict resolution paper ##
## this script will generate tables ##
## and conduct analysis ##
## 05.10.2013 -- df ##

## note, we were able to add a few dozen bills to our data since publication ##
## we have included these in the replication data, therefore the figures will not match exactly ##
## however there are no substantive differences in the model results ##

## clean up ##

	rm(list = ls(all = TRUE))

## load packages  ##

	library(stringr)
	library(Hmisc)
	library(mnormt)
	library(lubridate)
	library(lme4)
	library(survival)
	library(apsrtable)
	library(car)
	library(pscl)

## load data  ##

	data <- read.table('/Users/david/Dropbox/accepted/prq2013davidThomasOli/data/replication.txt', sep='|', header=T)

## reduce data to bills seen for the first time by br ##
## this reduces bias from "learning effects" ##

## check sample size to track reductions at each step ##

	cat("orginal data size: n = ",length(data[,1]))

## reduce to only government introduced bills ##
## or coalition member bills ##

	data <-  data[data$government == 1, ] 
	cat("n = ",length(data[,1]))

## cull sample of international accords ##

	droplist <- c(grep("x", data$GestaNr), grep("X", data$GestaNr))
	data <- data[-droplist, ]
	cat("n = ",length(data[,1]))


## how many bills, how many conciliation references ##

	{
	cat("There are ", length(data$conciliation), " total bills. \n", sep="")
	cat("Of these, ", length(data$conciliation[data$conciliation==1]), " went to conciliation. \n", sep="")
	}


## *********** ##
## values for table 1 ##
## *********** ##

## first create buckets for monthly break down ##
## give each bill a 'bucket' for histogram ##
## 'bucket' is simply the number of months ##
## elapsed since jan 1974 ##

	data$bucket <- c()

	for (i in 1:length(data$initiate)){
		data$bucket[i] <- ((data$iYear[i]-1974)*12) + data$iMonth[i]
	}
		
## mandatory gov ##
		length(data$mandatory[data$mandatory == 1 & data$brConc == 'gov'])
		length(data$mandatory[data$conciliation == 1 & data$mandatory == 1 & data$brConc == 'gov'])/length(data$mandatory[data$mandatory == 1 & data$brConc == 'gov'])

	## mandatory divided opp ##
		length(data$mandatory[data$mandatory == 1 & data$brConc == 'neu'])
		length(data$mandatory[data$conciliation == 1 & data$mandatory == 1 & data$brConc == 'neu'])/length(data$mandatory[data$mandatory == 1 & data$brConc == 'neu'])

	## mandatory unified opp ##
		length(data$mandatory[data$mandatory == 1 & data$brConc == 'opp'])
		length(data$mandatory[data$conciliation == 1 & data$mandatory == 1 & data$brConc == 'opp'])/length(data$mandatory[data$mandatory == 1 & data$brConc == 'opp'])

	## non-mandatory gov ##
		length(data$mandatory[data$mandatory == 0 & data$brConc == 'gov'])
		length(data$mandatory[data$conciliation == 1 & data$mandatory == 0 & data$brConc == 'gov'])/length(data$mandatory[data$mandatory == 0 & data$brConc == 'gov'])

	## non-mandatory divided opp ##
		length(data$mandatory[data$mandatory == 0 & data$brConc == 'neu'])
		length(data$mandatory[data$conciliation == 1 & data$mandatory == 0 & data$brConc == 'neu'])/length(data$mandatory[data$mandatory == 0 & data$brConc == 'neu'])

	## non-mandatory unified opp ##
		length(data$mandatory[data$mandatory == 0 & data$brConc == 'opp'])
		length(data$mandatory[data$conciliation == 1 & data$mandatory == 0 & data$brConc == 'opp'])/length(data$mandatory[data$mandatory == 0 & data$brConc == 'opp'])

## t test ##

	mean(data$conciliation[data$brConc!='gov'])/mean(data$conciliation[data$brConc=='gov'])
	
	t.test(data$conciliation[data$brConc=='gov'], data$conciliation[data$brConc!='gov'])

## created folded committee composition measure ##

	data$foldNeg <- c()
	
	for (i in 1:length(data$concConc)){
		if(data$conConcDisp[i]=='hung' | data$conConcDisp[i]=='opp' | data$conConcDisp[i]=='govopp' | data$conConcDisp[i]=='hungopp' | data$conConcDisp[i]=='govhung'){
			data$foldNeg[i] <- 1
		}else{
			data$foldNeg[i] <- 0
		}
	}

## we have to drop two disposition sets due to insufficient variation ##
## this does not change the substance of the results ##
## it merely eliminates two meaningless coefficients ##

	data <- data[data$brConcDisp != 'oppgov' & data$brConcDisp != 'oppneu', ]
	length(data$brIni)
	data$govneu <- 0
	data$govopp <- 0
	data$neu <- 0
	data$neuopp <- 0
	data$opp <- 0

	data$govneu[data$brConcDisp == 'govneu'] <- 1
	data$govopp[data$brConcDisp == 'govopp'] <- 1
	data$neu[data$brConcDisp == 'neu'] <- 1
	data$neuopp[data$brConcDisp == 'neuopp'] <- 1
	data$opp[data$brConcDisp == 'opp'] <- 1
	
	table(data$brConc, data$concConc)

## logit model for Table 2 ##

	tab4 <-glmer(conciliation ~ mandatory + brConcDisp + foldNeg + (1 | Legislaturperiode), family = binomial(link = 'logit'), data = data)
	summary(tab4)

## check for 2006 reform differences ##

	ref2006 <- glmer(conciliation ~ mandatory + brConcDisp + foldNeg + (1 | Legislaturperiode), family = binomial(link = 'logit'), data = data[data$jInitiate < 13122, ])
	summary(ref2006)

## no substantive cahnge ##
######################

## fit check ##

	empty <- glm(conciliation ~ 1, family = binomial(link = 'logit'), data = data)
	summary(empty)

## deviance ##
	dev <- deviance(tab4) - deviance(empty)
	pchisq(dev, 1, lower.tail = TRUE)

## log likelihood ##
	chi <- 2*(logLik(tab4) - logLik(empty))[1]
	pchisq(chi, 9, lower.tail = FALSE)
	
## mcfadden ##

	sqrt(1-(logLik(tab4)/logLik(empty)))

## probabilities for Table 3 ##

	simCoef <- rmnorm(n = 10000, mean = fixef(tab4), varcov = matrix(vcov(tab4), ncol = 8))
	govMan <- exp(simCoef[, 1] + simCoef[, 2])
	govNm <- exp(simCoef[, 1])
	mixMan <- exp(simCoef[, 1] + simCoef[, 2] + simCoef[, 5])
	mixNm <- exp(simCoef[, 1] + simCoef[, 5])
	govMixMan <- exp(simCoef[, 1] + simCoef[, 2] + simCoef[, 3])
	govMixNm <- exp(simCoef[, 1] + simCoef[, 3])
	hosMan <- exp(simCoef[, 1] + simCoef[, 2] + simCoef[, 7])
	hosNm <- exp(simCoef[, 1] + simCoef[, 7])
	mixHosMan <- exp(simCoef[, 1] + simCoef[, 2] + simCoef[, 6])
	mixHosNm <- exp(simCoef[, 1] + simCoef[, 6])

## put it out ##

	quantile(govMan, .5)
	c(quantile(govMan, .025), quantile(govMan, .975))
	quantile(govNm, .5)
	c(quantile(govNm, .025), quantile(govNm, .975))
	quantile(mixMan, .5)
	c(quantile(mixMan, .025), quantile(mixMan, .975))
	quantile(mixNm, .5)
	c(quantile(mixNm, .025), quantile(mixNm, .975))
	quantile(govMixMan, .5)
	c(quantile(govMixMan, .025), quantile(govMixMan, .975))
	quantile(govMixNm, .5)
	c(quantile(govMixNm, .025), quantile(govMixNm, .975))
	quantile(hosMan, .5)
	c(quantile(hosMan, .025), quantile(hosMan, .975))
	quantile(hosNm, .5)
	c(quantile(hosNm, .025), quantile(hosNm, .975))
	quantile(mixHosMan, .5)
	c(quantile(mixHosMan, .025), quantile(mixHosMan, .975))
	quantile(mixHosNm, .5)
	c(quantile(mixHosNm, .025), quantile(mixHosNm, .975))

## `accomodative conciliation committee' test ##

	table(data$negConPivot)
	table(data$negPivot)
	data$wild <- data$posConPivot - data$posPivot
	
	t.test(data$brCon[data$wild == 0], data$brCon[data$wild == 1])

## hierarchical model ##

	hier <-glmer(brCon ~ mandatory + wild + (1 | conConcDisp) + (1 | brConcDisp) + (1 | Legislaturperiode), family = binomial(link = 'logit'), data = data)
	summary(hier)
